# chirp
# Tim Menzies
# here

# todo: compute median once at end of clustering
# find border points of a cluster- whare at least one neighbor is
#     not in the cluster
# radiate out from the cluster, one grid at a time, break ties by
#     size of nearest cluster. which actually means sort clusters
#     by size, the from biggest to smallest, extend border by 1

@include "lib.awk"
@include "dist.awk"
@include "readcsv.awk" 

function _chirp(  o,_Chirp) {
  if (args("-d,data/nasa93.csv,-m,3,--norm,1,-seed,1,-cells,128",o)) {
    m     = o["-m"]
    normp = o["--norm"] 
    cells = o["-cells"]
    chirp(o["-d"],o["-seed"],_Chirp) } 
}
function chirp(file,s,_Chirp,_Rows,   one) {
  seed(s)
  readcsv("cat " file,_Rows)
  if (cells > length(d)) 
    cells = length(d)
  div = round(log(cells)/log(m))
  print "div1", div
  div = length(d) / div
  print "div2", div
  #div = 3
  one = corner1(any(all),_Chirp,_Rows)
  corners[1] = one
  sides(1,one,_Chirp,_Rows) 
  grid(_Chirp,_Rows)
}
function reset_Grid(_Grid) {
  a(id); a(box); a(ats); a(n); a(s); a(s2); 
  a(cluster); a(owner); a(descending); a(neighbor)
}
function new_Grid_Id(_Grid,     another) {
  another     = length(id)+1
  id[another] = another
  return another
}
function grid(_Chirp, _Rows, _Grid) {
  reset_Grid(_Grid)
  useds(_Chirp,_Grid)
  initClusters(_Chirp,_Rows,_Grid)
  sortClusters(_Grid)
  look4Neighbors(_Chirp,_Grid)  
  gridCluster(_Chirp,_Grid,_Rows)
}

function initClusters(_Chirp,_Rows,_Grid, at,cl) {
  for(at in ats) {
    cluster[++cl][1] = at
    addPoints2cluster(box[ats[at]],cl,
		      _Chirp,_Rows,_Grid) }
}
function sortClusters(_Grid) {
  asort(descending,descending,"xsortDown") 
  showClusters(_Grid)
}
function gridCluster(_Chirp,_Grid,_Rows,  
		     i,unused,max,again,holds) {
  for(i in descending)
    unused[ ++max ] = 1 
  do {
    again = 0
    print 1
    for(i=1;i<=max;i++) 
      if (unused[i]) {
	unused[i] = 0  # start from here and work right
	again = gc1(i,max,unused,holds,
		    _Chirp,_Grid,_Rows) }
  } while (again)
}
# to do : back out if spread increases
function gc1(i1,max,unused,holds,_Chirp,_Grid,_Rows,
	     cl1,at1,i2,cl2,at2,again,near,this,point,points,klass) {
  a(points)
  klass = klass1(gp)
  cl1 = descending[i1]["="]
  at1 = cluster[cl1][1]
  print "\n[" (holds[cl1][1] = at1) "]"
  spread[cl1] = morePoints(at1,points,_Grid,_Rows)
  med[   cl1] = medianOfPoints(points)
  owner[at1] = cl1
  unused[i1] = 0
  for(this in neighbor[at1]) near[this]
  for(i2=i1;i2<=max;i2++) # work from here down 
    if (unused[i2]) {
      cl2 = descending[i2]["="]
      at2 = cluster[cl2][1]
      if (at2 in near)
	  if (canJoin(cl1,cl2,_Grid)) {
	    for(this in neighbor[at2]) near[this]
	    again = 1
	    unused[i2] = 0
	    push(at2, holds[cl1])
	    spread[cl1] = morePoints(at2,points,_Grid,_Rows)
	    med[   cl1] = medianOfPoints(points)
	    owner[at2] = cl2
	    n[cl1] +=  n[cl2]
	    s[cl1] +=  s[cl2]
	    s2[cl1] += s2[cl2] 
	    print "adding " at2 " to " holds[cl1][1], 
	          " :n " n[cl1]" :mu "s[cl1]/n[cl1]" :sd " sd(n[cl1],s[cl1],s2[cl1]),
	           " :med " med[cl1] " :spread " spread[cl1]
	  }}
  return again
}
function medianOfPoints(lst) {
  return lst[  int(length(lst)*0.5) ]
}
function morePoints(at,points,_Grid,_Rows,    point,klass,max) {
  klass = klass1(gp)
  for(point in box[ats[at]]) { 
    point = box[ats[at]][point]
    push(d[point][klass],points)
  }
  max = asort(points)
  return points[ int(max*0.75) ] - points[ int(max*0.25) ]
}
function useds(_Chirp,_Grid,  
	     point,at,t) {
  for(point in xndex)  { # for all the points
    at = used(point,_Chirp,_Grid) # ensure we know this at
    if (at in ats)
      t = ats[at] 
    else
      t = new_Grid_Id(_Grid)
    ats[at] = t 
    push2(t,point,box)   # remember we used them
}}
function used(point,_Chirp,_Grid,   at,sep,pole) {
  # as a side effect of locating a point,
  # we also remember its n-dimenaional location
  # (where it is "at"
  for(pole=1;pole<=m;pole++) {
    at = at sep xndex[point][pole]
    sep = ","
  }
  for(pole in xndex[point]) 
    axis[at][pole] = xndex[point][pole]
  return at
}
function addPoints2cluster(lst,cl,_Chirp,_Rows,_Grid,  
			   i,point,dep,klass) {
  klass = klass1(gp)
  for(i in lst) {
    point   = lst[i]
    dep     = d[point][klass]
    n[cl]  += 1
    s[cl]  += dep
    s2[cl] += dep*dep
  }
  descending[cl]["x"] = n[cl]
  descending[cl]["="] = cl
}
function showClusters(_Grid,   i,j,max) {
  max = length(descending)
  for(i=1;i<=max;i++)  {
    j = descending[i]["="]
    print(" :birth ["cluster[j][1]"] :id " j" :n "n[j] \
          " :mean " s[j]/n[j] " :sd " sd(n[j],s[j],s2[j]))
    }
}

function canJoin(i,j,_Grid,  old,new) {
  old = sd(n[i],s[i],s2[i])
  new = sd(n[i] + n[j],
	   s[i] + s[j],
 	   s2[i]+s2[j])
  return new <= old  
}
function lookThisWay(way,_Chirp,_Grid,
		  at1,at2,sep,pole,old,new) {
  for(at1 in ats)  {
    at2=sep=""
    for(pole=1;pole<=m;pole++) {
      old = axis[at1][pole]
      new = old + way[pole]
      if (new <    0) new = 0
      if (new >= div) new = div - 1
      at2 = at2 sep new
      sep = ","
    } 
    if (at2 != at1) 
      if (at2 in ats)  { # if it is a used grid 
	neighbor[at1][at2] = 1
	neighbor[at2][at1] = 1}}}

function look4Neighbors(_Chirp,_Grid,
		   lo,hi,i,way,top) {
  hi = 1 ; lo = -1
  for(i=1;i<=m;i++) 
    way[i] = lo
  while (1) {
    lookThisWay(way,_Chirp,_Grid)
    top=1
    for(i=1;i<=m;i++) 
      top = top && way[i] == hi
    if (top) 
      break
    way[m]++             # set the last way
    for(i=m;i>=1;i--)    # ripple left any overflow
      if (way[i] > hi) {
	way[i] = lo
	way[i-1]++
	if (way[i-1] <= hi)
	  break }}}


function corner1(corner0,_Chirp,_Rows,   \
                 tmp,i,most,out) {
  for(i in all)  {
    tmp = Dist(corner0, all[i],_Chirp,_Rows)
    if (tmp > most) { 
      most=tmp; out = all[i]}}
  return out
}
function sides(p,one,_Chirp,_Rows,    two) {
  if (two = furthest(_Chirp,_Rows)) {
    print p, pole(p,one,two,_Chirp,_Rows)
    corners[++p] = two    
    if(p <= m) 
      sides(p,two,_Chirp,_Rows)}
}
function furthest(_Chirp,_Rows,  \
                  i,one,two,tmp,most,out) {
  for(i in all) {
    one = all[i]   
    if (has(corners,one)) continue
    tmp = 0 
    for(two in corners)
      tmp += Dist(one,corners[two], _Chirp,_Rows)
    if (tmp > most) { out = one; most = tmp}
  }
  return out 
}
function pole(p,one,two,_Chirp,_Rows,\
              most,least,i,three,aa,bb,cc,xx,yy,what,where) {
  poles[p] = p
  l[p] = one
  r[p] = two
  c[p] = cc = Dist(l[p],r[p],_Chirp,_Rows)
  most = -1 * (least = 10**32)
  for(i in all) {
    three = all[i]
    if ( three == l[p] ) {
      aa = xx = 0; bb = cc
    } else {
      if (three == r[p] ) {
        aa = xx = cc; bb = 0
      } else {
        aa = Dist(three, l[p], _Chirp, _Rows)
        bb = Dist(three, r[p], _Chirp, _Rows)
        xx = (aa^2 + cc^2 - bb^2)/(2*cc + 0.0000001)
      }}
    x[p][i]  = xx
    yy = sqrt(abs(aa^2 - xx^2)) 
    y[p][i]  = yy
    sorted[p][i]["="] = three 
    sorted[p][i]["x"] = xx
    if (xx < least) { lhs[p] = three; least = xx} 
    if (xx > most ) { rhs[p] = three; most  = xx}
    err[p]  += yy/length(all)
  }
  asort(sorted[p],sorted[p],"xsort")  
  for(i in sorted[p]) {
    what = sorted[p][i]["="] 
    where = int(i/div) 
    if (where > m) where--
    xndex[what][p] = where
  }
  return err[p]
}
function print_Chirp(_Chirp,   pmax,p,xmax,i) {
  pmax = length(l)
  print pmax
  for(p=1;p<=pmax;p++)  {
    print ""
    xmax = length(sorted[p])
    print p" :left "lhs[p]
    for(i=1;i<=xmax;i++)
      print p,i" :j "sorted[p][i]["="]" :x "sorted[p][i]["x"]
    print p" : right "rhs[p]
}}
function Dist(i,j,_Chirp,_Rows,   out) {
  if ( (i,j) in cache ) 
    return cache[i,j]
  out = dist(i,j,_Rows,normp)
  cache[i,j] = cache[j,i] = out
  return out
}
